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Abstract 


In dynamic fracture problems role of material inertia becomes significant.Tlie for- 
mulation and analytical solutions of dynamic fracture problems are of very complex 
nature. Experimental and numerical methods like finite element method are used, 
extensively for solving dynamic fracture problems. Existing methods for simulating 
crack phenomenon in dynamic fracture problems such as, moving mesh procedures 
or force release models, have certain drawbacks such as, they are either computa- 
tionally very intensive or give oscUlations in the solution. 

In this work, a model is proposed for simulating dynamic crack propaga- 
tion phenomenon by finite element method. Gradual propagation of crack from one 
node to the next node is simulated by adding a 1-D elastic element at the crack tip 
node, whose stiffness value is reduced gradually to zero as the crack propagates to 
the end of next element. 

The finite element model is first validated for static case and then with 
hypothetical data for different crack speeds. Results give smooth variation of en- 
ergy release rate curve for various crack velocities. This method is also applied to 
determine the dynamic Gi for the experimental data. Results show that the energy 
release rate gives a solution. 
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Chapter 1 

INTRODUCTION 


1.1 Introduction 

The dynamic fracture phenomenon can be characterized by various dynamic states 
of crack tip. The dynaunic states of crack tip are induced by impaict loading applied 
to the cracked soUds, or by fast motions of the crack tip itself. Dynamic loads may 
be created by moving vehicles, wind gusts, shocks or blasts, unbalanced machines, 
wave impacts, seismic disturbances etc. For solving these problems often the time 
variation is neglected and they are treated as static or quasi-static problems. Struc- 
tural members that are subjected to cyclic loading such as vibrations, microscopic 
imperfections in the material give rise to microscopic cracks, which in turn grow into 
cracks of significant size over passage of time. At this stage further growth of crack 
in a stable or unstable manner depends not only on the material but also on the 
nature of loading, i.e static or dynamic. For many dynamic problems it is impossible 
to find the closed form solutions, especially if they involve the aspect of fracture. A 
number of analytical solutions to problems of dynamic fracture shedding important 
light on the basic phenomenon have been obtained in the past two decades. These 
solutions are however limited to simple cases of loadings and of unbounded plane 
bodies. Usuallythe interactions of the stress waves emanating from the crack tip 
and of those reflected from the boundaries of a finite body make the analytical solu- 
tions invalid. Thus, it often becomes mandatory to analyze dynamic crack problems 
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in finite solids using numerical method. 

The subject of the dynamic fracture can be broadly described as that of 
mechamcs of solids, wherein the effect of material inertia and stress waves on crack 
play significant roles. One method of classifying fracture problems is as follows 
(Nishoka and Atluri, [1986]): 

1. Solids containing stationary cracks subjected to dynamic loading 

2. Solids containing dynamically propagating cracks under quasi-static loading 

3. Solids containing dynamically propagating cracks under dynamically loading 

Dynamic fracture mechanics concerns with (i) the onset of crack growth 
under dynamic loading (ii) dynamic crack propagation in stressed solids. The con- 
ventional linear elasto-plastic (or) quasi-static fracture mechanics applies ordy up 
to the end of stable growth under quasi-static loadings and assumes that the on- 
set of unstable crack propagation renders the structure useless. Yet the prevention 
of crack growth initiation itself may be too conservative or too costly an objective 
for some structural designs, and also catastrophic failures caused by unstable crack 
growth are obviously intolerable. In these cases, the assurance of crack arrest, as 
a second line of defense, is essential. This concept has been investigated for design 
methodologies assuring the integrity of nuclear pressure vessels under thermal shock 
conditions, LNG ship hulls, and gas transmission pipelines. 

The fundamental frame work of the subject of dynamic fracture mechanics 
relies primarily on the solutions for dynamic behaviour of solids containing cracks. 
These solutions are characterized by (i) stress integrations (ii) the need to account 
for kinetic energy in the global energy balance of the fracturing body and (iii) inertia 
effects of the material. The inertia effect may arise due to two reasons (i) rapidly 
applied load on a cracked solid and (ii) rapid crack propagation. In the first case 
the influence of the load is transferred to the crack by stress waves through the 
material. In the second case the material particles on the two crack faces displace 
each other, as the crack advances. The inertia effect, in the first case, is considered 
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significant when the tune taken to load the specimen to maximum value is small 

as compared to the time required for a characteristic stress wave of the material 

\ 

to travel a characteristic dimension of the body. In the second case, the inertia 
effect should be accounted for whenever the velocity is a significant fraction of the 
characteristic velocity(e.g. Rayleigh wave velocity) 

In the past two decades or so, a number of eineJytical solutions,which pro- 
vides a useful understanding of dyn ami p crack behaviour, have been obtained. These 
analytical solutions are however, limited to cases of simple loadings and unbounded 
plane bodies. Moreover the stress wave introductions, which play an important role 
in dynamic fracture mechanics usually render the analytical solutions intractable. 
Therefore, the use of numerical methods is often indispensable for the analysis of 
cracks in finite solids. 

1.2 Literature Survey 

Though considerable amount of work has been done to study the fracture phenom- 
ena through numerical methods, only large rectangular cantilever beam specimens 
were given special interest. Owen and Shantaram [1977] started the case of finite 
element method to study dynamic crack growth. They studied double cantilever 
beam specimen and pipeline problems under transient loading. The crack was ad- 
vanced from one node to the next node, when the stress at the gauss point nearest 
to the crack tip exceeded a certain value and then damping coefficient was made 
zero at the released node. The crack propagation history was simulated for the given 
loading conditions. 

Nishoka and Atluri [1982b] investigated the crack propagation and arrest 
in a high strength steel DCB specimen using moving singular dynamic finite ele- 
ment procedure. An edge crack in a rectangular DCB specimen was propagated by 
inserting a wedge. The results were compared with experimental data of caustics. 
In another work, Nishoka and Atluri [1982a] presented the results of generation and 
prediction studies of dynamic crack propagation in plane stress and plane strain 
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cases. The studies were conducted by using FEM, taking into account stress sin- 
gularity near the crack tip. The variation of dynamic stress intensity factor with 
time and the variation of dynamic fracture toughness with velocity were studied and 
compared with available experimental results. 

Nishoka and Atluri [1986] gave elaborate information on analysis of dy- 
namic fracture using FEM. The most common ways to deal with the crack tip region 
are to simulate crack growth through gradual release of element nodal forces or em- 
bedding a moving element in which the interpolation functions are determined by 

description the differential equations in time for the node point variables must be 
integrated. Because the dynamic fields associated with rapid crack growth are rich 
in high frequency content, small time steps are required for accuracy. Because of 
natural restriction to small time steps, the authors report that many times combina- 
tion of conditionally stable explicit time integration scheme and diagonalised mass 
are found to be accurate. 

Chiang [1990] presented a numerical procedure based on eigen functions 
to determine the dynamic stress intensity factor of crack moving at steady state 
under anti-plane strain condition. An edge problem and a radial crack problem 
were solved using traction and displacement boundary conditions separately. It was 
shown that the dynamic effect is relatively insignificant for low crack propagation 
speeds provided that specimen size is fairly large. 

Thesken and Gudmundson [1991] worked with an elasto-dynamic moving 
element formulation incorporating a variable order singular element to enhance the 
local crack tip description. Kennedy and Kim [1993] incorporated micropolar elas- 
ticity theory into a plane strain finite element formulation to analyze the dynamic 
response of the crack. Materials with strong micropolar properties were found to 
have significantly lower dynamic energy release rate than their classic material coun- 
terparts. Wang and WilUans [1994] investigated high speed crack growth in a thin 
double cantilever beam specimens using FEM. The cantilever end was loaded un- 
der an initial step displacement of 1 mm and then pulled with a constant velocity 


the continuum neax tip in the mesh, and J-Integral considerations. FoUowing ^acia l 
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and the crack was propagated at assumed speed by gradual node release technique. 
Large dynamic effects were observed because of wave reflections within finite speci- 
men size. The reflection and interaction of the waves from the free boundary were 
avoided by limiting the duration of study. Beissel, Johnson and Popelar [1998] pre- 
sented an algorithm which allows crack propagation in any direction but doesn’t 
require remeshing or the definition of new contact surfaces. This is achieved by 
tracking the path of the crack tip and failing the element crossed by the path such 
that they can no longer sustain tensile volumetric stress. The edges of these faded 
elements simulate crack faces that can sustain only compressive normal traction. Lin 
and Smith [1997] presented a method wherein crack growth is predicted in a step by 
step basis from Paris law using stress intensity factor calculated by using finite ele- 
ment method. The crack front is defined by a cubic spUne curve from a set of nodes. 
Both the one quarter-node crack opening and 3-D J-Integral method were used to 
calculate stress intensity factor. Automatic remeshing of finite element model to 
a new position that defines the new crack front enables the crack propagation to 
be followed. Lee, Hawong and Choi [1996] studied propagating crack problems of 
orthotropic material under the dynamic plane mode. Dynamic stress components 
and dynamic displacement components moimd the crcick tip of an orthotropic mate- 
rial under the dynamic load and the steady state in crack propagation were derived. 
When the crack propagation speed approaches zero, dynamic stress components and 
dynamic displacement components are identical to those of static state. The stress 
component values of the crack tip are greater when the fiber direction coincides with 
the direction of the stress component than when the fiber direction is normal to the 
stress component. Chandra and Krauthammer [1995] investigated effects of rapidly 
applied load on the J -integral and stress intensity factor (AT ) of a cracked solid. 
Variation of J has been studied using energy balance approach, whereas that of K 
has been dealt with elastodynamic considerations together with several simplifying 
assumptions and approximations. 

Christina and Christer [2001] presented a method for obtaining the com- 
plex stress intensitv factor for an interface crack in a bunaterial using mmimum 
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number of computations. A crack closure integral method has been used. Nageswara 
Rao and Acharya [1995] studied the mode I interlaminax fracture toughness, G/c 
of composites using a slender double cantilever beam specimen. Zhuang 8ind Guo 
[1999] addressed recent developments in the area of dynamic fraicture mechanics and 
the applications of analysis methods to the rapid crack propagation for gas pipelines. 
Criteria for crack initiation, propagation and arrest were discussed. 

1.3 Present Work 

Present work is a development of ‘StiflEness Release Model’ for dynamic crack prop- 
agation in DCB specimen by Si\n. Reddy [ 1997 ]. An additional one dimensional 
elastic element is attached at the crack tip node whose stiffness value is gradually 
reduced to zero as the crack advances to the end of next element. The thesis is 
organised as follows; 

Chapter 2 describes the basis of FEM formulation and existing propaga- 
tion model. 

Chapter 3 deals with the crack propagation element model. 

Chapter 4 presents the validation of code and results and 

Chapter 5 presents the conclusion of the present work and scope of future 

work. 
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Chapter 2 


FEM FORMULATION AND 
CRACK MODELLING 

2.1 FEM Formulation 

The finite element solution involves discretization of domain (Q) into suitable ele- 
ments and approximating the field variable interior to the element in terms of its 
nodal values using suitable shape function. 

Equations that govern the dynamic response of a system are derived by 
using principle of minimum potential energy or principle of virtual work. Principle 
of virtual work states that ‘If a set of stresses is in equilibrium at all points with a 
set of external tractions, then the sum of internal and external work done during 
any virtual displacement must be zero.’ 

In the following paragraphs, FEM equations are derived in brief using 
matrix notation. If {Pb} is a set of body forces acting over the domain, Cl and 
{F 5 } is a set of traction acting over the boundary (Fo-). Then, stresses and strains 
induced are {a} and {e} respectively. Then, 

Work done by body force = 

n 

Work done bv surface traction = /{dg}{Ps}d5 
Internal work done by stresses = 

Internal work done by inertial forces = J {5 g}'^ p{q} dV 
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where 


Internal work done by damping forces = S{3q}^kd{q}dV 

n 

{<^ 9 }, {(J<} : Virtual, displacements and corresponding strains 

{sill?} • Velocity and acceleration component 

{cr} : Stress component 

p ; Density of material 

kd : Damping coefficient 

{Pb} ■ Body forces 

{Ps} '■ Surface traction 

Tff : Surface traction boundary 


The virtual work principle i.e., 

sum of internal work = sum of external work. 

leads to, 

J{S€}^{a}dV-t- l{Sqyp{q}dV + J lSq]^ki{q}dV 

n n n 

= / {5qf{PB}dV + I {5qf{Ps}dS (2.1) 

n To- 

In F.E. method, the displacements within element are expressed in terms of 

its nodal values as 

{g}(‘) = (2.2) 
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where, 


[iV](®) - Shape function matrix 

Strain-displacements relation can be written as 

(2.3) 

where, 

{e} - Strain component 

- Strain displacement matrix 

Stress-strain relation and stress-nodal displacement is given as 

= (2.4) 


where, 

- Material property matrix 

The present analysis is made for an undamped system and neghgible body 
forces. The displacements, strains and stresses axe expressed in terms of nodal 
variables as given in equations 2.2, 2.3, 2.4. This on substitution in Eqn 2.1 gives. 


^ ( r \ 

E [ {g} + E / {g}- 

/=•' vji) J / 


Ne r 


where summation is over all elements, Ne- 
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Since {(5g} is to be arbitrary, satisfying kinematic conditions, the equation can be 
simplified as, 


W-> / <«) U.) 


(«) 
Ne 

= E 

e=l 




\ 


I {PsYMd^drj 

Vi'’ 


2 . 6 ) 


where, 

1 J] determinant of Jacobian matrix. 

which can be written as 

Ne Ne 

(2-7) 

e=l c=l c=l 

where 

[K](e) = / [D]^^^[B]^^^\J\d^dT ] , Elemental stiffness matrix 

n(«) 


Jlci^dr; , Elemental mass matrix 

n(-=) 


= j {PsY‘‘^J\d^d,ri, Elemental force vector 

r(<) 

The final set of assembled equations can be written as. 




( 2 . 8 ) 


where 

[A/] = ^[A/]^®\ Global mass matrix 
[A'] = Global stiflfiiess matrix 
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{F} — Global applied force vector 


2.1.1 Integration algorithm 

Equilibrium equation 2.8 can be solved either by time integration or by mode su- 
perposition of which former one is preferred in the wave propagation problems. In 
this integration scheme, there are many different methods, which can be classified 
as “Explicit” or “Implicit” which has specific advantages and disadvantages ( Bathe 
[1990] ). In the present work two different integration schemes are used 

1. Newmark’s method 

2. 9 - Method 

which are discussed in the following section. 

Newmark’s Method of Integration 

Newmark’s method is an extension of ‘linear acceleration’ method, in which, accel- 
eration is assumed to be linear within each time step, t to t + At. 

Thus, velocities and accelerations are expressed in terms of displacements as, 

= {«}■ + [(i- < 5 ){ 9 }‘ + ( 2 . 9 ) 


«} 


i-f- At 


= {,}• + {«}‘Ai + [(i - a) {(,•}< + <S{gr^' 


Af 


( 2 . 10 ) 


where a and /? are the parameters chosen smtably to have accuracy and stability. 
Usually, oc and p are taken as 0.25 and 0.5 respectively to get unconditional stability. 
Solving equations 2.9 and 2.10, we get expressions for and in terms 

unknown displacements These are then substituted in the the following 

equation to solve for 

[M]{qY^^^ + = {FY^^' ( 2 . 11 ) 
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The above method can be summarized as follows. ( Bathe [1990] ) 

1. Form global stiffiaess matrix [K], and mass matrix [M]. 

2. Initialize {9}°, {5}°, {5}°. 

3. Select time step At and then calculate integration constants 


ao 


1 

OiAt^ ’ 


ai 


6 

olAV 




1 





Ofi = At(l - 5); 


(I7 = 5 At 


Form eflFective stiffness, [K'] = [jK] + oo[M] 


4. For each step: 


(a) Calculate effective load at t + At 

{FT^‘ = + [M] + a,{5}‘ + a3{5}‘) (2.12) 

(b) Obtain the displacements at time t + At by = {F'} 

(c) Calculate acceleration and velocity at time t + At 

{«}'+" = ({«}‘^“ - {«}■) - (2.13) 

= {q}‘ + a6{?}‘ + (2.14) 
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0-Method 


HofF and Pahl [ 1986 ] presented an unconditionally stable second-order-accurate 6 
method. This method is found to be more stable compared to Newmark’s method 
which can be summarized in brief as follows. 

1 . Choose parameters di within 0.95 to 1.0 and 6q as 1.0 

2 . Form global stiffness matrix [K], and mass matrix [M]. 

3. Initialize { 5 }°, {q}^. 

4. Select time step At and calculate integration constants 

4.0 X Of (1.5 — 0i) X 4.0 X Of 

= At 

For effective stiffness, [K'^ = [K] -I- ao[M] 

5. For each time step : 

(a) Calculate effective load at t + At 

= {F}**-*' + [M] (a,({g} + Af{g} + 0.5 x At^s})' - {«"}) 

+(1.0-«i)At[i(n{g}‘-“ (2.15) 

(b) Solve for displacements at time t + At by 

(c) Calculate acceleration and velocity at time 1 4- At 

{,}•+“ = oc ({?}“+“ - {g}‘) - ooAt X {«}• + (l.O - 20?) {«}• (2.16) 

= {g}‘ + At («1 - 0.5) {ft' + (1.5 - 0i) At{ft' + At (2.17) 
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^-Method offers an unconditionaUy stable solution. WMle determining resultant 
force from the displacement boundary condition, Nevnnaxk’s method was observed 
to give imstable oscillations . On the other hand, d~ method was found to give stable 
solution under similar conditions. 

2.2 Methods for Simulation of Dynamic Crack 
Propagation by FEM 

To simulate a crack propagation in soUds two different concepts of finite element 
modelling are in use i.e., stationary mesh procedure and moving mesh procedure. 
Out of these two, moving mesh procedure, as the name implies involve change of 
mesh in each step as the crack propagates. This is computationally very intensive. 
Stationary mesh procedure do not require changes in mesh, and require only change 
in the boundary conditions. 


2.2.1 Stationary Mesh Procedure: 

In the simple stationary mesh procedure of modelling linear elastodynamic crack 
propagation, the nodes ahead of the crack tip are spaced at VcAt {Vc being crack 
velocity) and crack propagation is simulated by releasing one node at a time. How- 
ever, if a simple node release technique is used, release of constraint on the displace- 
ment (at the preceding crack tip) one in one time step amounts to discontinuous 
crack propagation and induces spurious highly oscUlatary finite element solution. To 
overcome these inaccuracies, several algorithms have been suggested in the hterature 
to release the node gradually over few time steps. 

Force Release Model 

Suppose that the actual tip is located at ‘C’ in between the finite element nodes 
B and D as shown in Fig 2.1. The length segments BC and BD are b and d re- 
spectively. The holding back force, F, at node B is gradually reduced to zero over 
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a number of time steps as the crack tip reaches to the node D. Various schemes 
available to decrease the force to zero as follows: 


1. Malluck and King [1978] suggested the release rate based on constant stress 
intensity factor. 


i-t' 


( 2 . 18 ) 


where Fo is the original reaction force when the crack tip was located at node 
B. 


2 . 


Rydhom et al [1978] suggested the release rate based on constant energy release 
rate. 



( 2 . 19 ) 


3. Kobayasi et al [1978] suggested the linear rate based on no physical argument 
other than the pure intuition. 


F 

Fo 



( 2 . 20 ) 


In order to have more gradual and smooth propagation of crack Kishore, 
Kumar and Verma [1993] used a modified method. Holding back force at the crack 
tip B is linearly decreased more gradually to zero when the crack tip reaches end of 
the next element. Thus when the crack tip goes beyond node B then 


= h_ 

Pub V d + d2 


(2.21) 


where is the reection 3 Lt node B, when the node was closed, b is the cra.ck 
extension and d^di are the elements length as shown in Fig 2.1. And when the crack 
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moves beyond node D to point Di 


_ f _ d d- 

Fhb V d + dj j 


( 2 . 22 ) 


i^ = fi 

Fhd V di + ^2 y 


(2.23) 


where Fhd the reaction at node D, when the node was closed, b, 62 the crack 
extension and d, di, d 2 are the elemental lengths as shown in Fig 2 . 1 . 


In force release model, time increment At is nsucilly taken such that it 
takes 15 to 20 iterations for the crack tip to cross one element. The amount of 
force necessary to keep the nodes together is determined cind a factor is used to 
proportionately decrease the force as the crack tip advances. Thus, at subsequent 
times, though the problem is highly dynamic, the current calculations were based 
on the force calculated to keep the nodes together 15 to 20 time steps before. This 
will cause discrepancy as the crack advance further from one element to the next 
one causing large oscillations in the solution. This discrepancy should be overcome 
for obtaining better solution. A new model is proposed to overcome the drawbacks 
of force release model which instead of dynamic force at the crack tip, takes into 
account the stiflhiess and mass of the element which is described in Chapter 3. 
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d 

di 

dl 





Figure 2.1: Crack opening scheme in force release model 
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2.3 Energy Releaise Rate and Its Determination 

In fracture mechanics a crack can be characterized by four parameters 

1. Energy release rate, Gi,Gn, Gm is energy based which, can be applied to 
brittle and less ductile material. 

2. Stress intensity factor (/C/, Kij,Kijj) is based on stress. This parameter is also 
applied to brittle and less ductile material. 

3. J-Integral (J) which has been developed to deal with ductile material which 
can also be applied to brittle material as well. 

4. Crack tip opening displacement (CTOD) is displacement based and is devel- 
oped for ductile material. 

From the computational point of view, stress intensity factor require spe- 
cial element so as to simulate stress singularity at the crack tip, whereas J-Integral 
and energy release rate do not require any such modelling of stress singularity. In 
energy release rate approach, energy is found out in entire domain and in J-integral, 
evaluation is done along a path far away from crack tip. Thus, both these approaches 
avoid analysis close to the crack tip and don’t need any modelling of stress singu- 
larity. 


Griffth realized that a crack in a body will not extend unless energy' is 
released in the process to overcome the energy needs of forming two new surfaces; 
one below and one above the crack plane, which is called surface energy of body- So, 
crack in a body will advance when the reduction in the total energy of the body is 
equals or more than the surface energy of the new surface thereby created. Most of 
the energy release, as the crack advances, comes from the parts of the plate which 

are adjacent to the cracked surface. 

The two important parameters need to be considered are, 


18 



1. How much energy is released when a crack advances, and 

2. Minimum energy require for crack to advance in forming two new surfaces. 


First parameter is measured in terms of energy release rate denoted by G. 
Which in general as a function of crack size. Energy release rate is the amount of 
energy release per umt increase in area during crack growth. This energy is supplied 
by the elastic energy in the body and by the loading system. 

The energy requirement for a crack to grow per unit area extension is .called as 
crack resistance and is usually denoted by symbol R. Crack resistance is the sum of 
energy required to (1) form new surface and (2) cause an elastic deformation. Both 
the available energy release rate and crack resistance are important to study the 
possibility of crack becoming critical. Obviously, when the available energy release 
rate far exceeds the crack resistance, the crack starts to grow at high speed. 

In the present study ‘energy release rate’ is adopted as it is a more com- 
prehensive concept. As the crack advances, 

1. StiflEness of the component decreases. 

2. Strain energy in the component either increases or decreases. 

3. Work is done on the component if there is an application of external load. 

4. Energy is being consumed to create two new surfaces. 


For quasi-static growth of crack by crack area AA, the incremental work 
by external force he A change in strain energy he A U, end the available en- 
ergy, GAA, satisfy energy balance equation. 

GAA = AW^-AU (2.24) 
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(2.25) 




G = - 


dU 

dA 


where 11 is the potential energy. 

For a plate of uniform thickness, dA = Bda 
where da is the incremental crack length 

B is the thickness of the plate 


G 


1 dU 
B da 


(2.26) 


(2.27) 


For dynamic crack propagation problem, the kinetic energy. T, of the body 
should be taken into consideration. Dynamic energy release rate, Go, is different 
from quasi-static case, as some energy may be consumed to impart kinetic energy 
to the cracked portion of the body and to generate stress waves. 

For dynamic case, energy balance becomes 

GAA = AWe^-AU-AT (2.28) 

where A Tis the increment in kinetic energy in the body. 

For constant velocity crack propagation, 

G = 4:^(1F^- r) (2.29) 

jd da 

For crack moving with constant velocity u 

Q ^ ^ (2.30) 

B V 
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In the present study, Energy release rate is found out using equation 2.29. 

In the Finite Element analysis different energy terms are found out as follows: 



(2.31) 


(2.32) 


(2.33) 


where 

Wa:t is the external work done 

U is the strain energy in the component 

Tis the kinetic energy in the component 

{F} is the external force vector 

{q} is the displacement vector 

{q}is the velocity vector 

[K] is the global stiffness matrix 

[M] is the global mass matrix 

2.4 Closure 

In this Chapter, formulation Of FEM equation is described along with 
integration algorithms for solving equiUbrium equation governing linear dynamic 
response of a system of finite element. Various methods of simulating dynamic 
crack propagation with stationary mesh are described. This Chapter also describes 
some general concepts in energy release rate which will be used in the next Chapters 
to construct a propagation element. 
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Chapter 3 

STIFFNESS RELEASE MODEL 


3.1 Introduction 

When the force release models were appUed to high speed crack propagation there 
were some drawbacks: 

1. The time increment, At is taken such that crack tip will take 15 to 20 time 
steps to cross one element. The amount of force necessary to keep the nodes is 
determined and is proportionately decreased as the crack advances. Though 
the problem is highly dynamic in nature, the current calculations were bcised 
on the force calculated to keep the nodes together at 15 to 20 time steps before. 
This will give rise to high oscillations as the crack advances from one element 
to next element. 

2. Further, to model the case of stopping of the crack propagation due to less 
available energy release rate, and starting again at a latter time, it would not 
be known how much holding force is to be applied. So, this models are not 
suitable to solve inverse problems. 

In stiffness release models above shortcomings can be overcome. 
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3.2 Stiffness Release Model 


In these stiffness release models, one-dimensional elastic element is added at the 
crack tip node, whose stiffness is reduced gradually as the crack advances. Spring 
stiffness is reduced gradually to zero as the crack reaches to the nest node ( Siva 
Reddy [1997]). 

Fig. 3.1 (a) shows sy mm etric part of a deformed plate in which crack is up to the node 
B. Fig. 3.1(c) shows the configuration of the plate when the crack has completely 
propagated to the next node C. Fig. 3.1(6) shows the modelling of plate with 
partially propagated crack (between B and C ) by the addition of one-dimensional 
linear spring element whose stiflEness is Ks- 

In brief, the development of the model consists of the following steps. 

Let Ko be the stiffness at node B of original plate without additional spring (i.e 
force required to cause unit displacement of node B) and iio be the displacement of 
the node when the crack tip reaches to the next node. 

Ks can be found out from the equilibrium equation as. 


Kq{uo -u) = KsU 

(3,1) 


(3.2) 


Now to express fCs as a function of crack size, it is necessary to know u as a 
function of crack size. For this purpose energy release rate is used. -A.t any stage, 
the crack can be extended by a small amount by reducing Ks by an amount AKs. 
This reduction in the spring causes a decrease in the total energy of the plate. As 
this is the only energy' loss in the system, this can be used while calculating energy 
release rate. 

The energy in the spring cm be wntten as 

E = \Ksn‘‘ 

= ( 3 - 3 ) 

2 V u J 
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The energy release rate is written as, 


G = 


AE 

AA 


where AA is increment in crack area. 

Let a non-dimensionalised parameter for crack length, ‘o’ be defined 
a = crack length (within the element) / d 
where d is the distance between two nodes. 

Therefore, the crack extension area, AA can be written as, 

AA = B dAa 

where B is the thickness of plate. 

For static case, energy release rate can be written as, 


G = 




1 du 
B dda 


Assuming linear variation of G in a, G can be written as 


G = Ki + K2a 


where Ki and K 2 are constants to be determined. 
Equating G from Eq.(3.5) and Eq.(3.6) 


„ „ IKqUo ( du\ 

Ki + K2a 2 5 ^ 


{^Ki ”1“ K 2 (x) B d doi — 2 ^0 ^0 du 

Integrating on both sides 

Bd 


K,a + ^ + C, 


= -■ r Ao UqU 


where Co is integration constant 


u = — 


Kq Uq 


-Bd 


Kl CL + 


Koa? 


+ Co 


By applying initial condition u = 0 wheno — 0, Co — 0 



using other end condition u = uo when a = 1.0, 


Uq = — 


Kq Uq 


Bd 




Eq. (3. 10) and Eq.(3.11) gives 


(3.11) 


u 

Uo 


On simplification it reduces to 


2 iiTi g + JCa 
2 Ki + K 2 


(3.12) 


— — Cl + C7(a — (z^) = f 

Uq 

where 

c = -f ) 

\2Ki + K2J 

Eq.(3.2) finally can be written as 


Ks = Kq 



(3.13) 


(3.14) 


(3.15) 


where / = u/uq as given by Eq.(3.13) 

Static stiffness Kq can be determined for a given specimen under static condition. 
In the above model stiffness was modified and mass of the element was maintained 
constant. Deshmukh [2000] proposed a model in which as the stiffness of the addi- 
tional elastic member at the crack tip decreases, the effective mass of the element 
which contain crack tip is increased as per modified shape function as the crack 
advances. 

These models overcomes the shortcomings of force release models, but when the 
crack moves from one element to next element there are oscillations in energy release 
rate curv'e, and it is not smooth when crack moves from one element to next one as 
expected. 
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(a) 














_A 

B 

C 


^ r 




^“0 ^ 

(c) 


Figure 3.1: Crack opening scheme in stifness release model 
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3.3 Present Stiffness Release Model 

In the present model, the additional spring stiffness Ks, is decreased to zero over 
two elements instead of one element, i.e., spring stiffness was gradually reduced to 
zero as the crack reaches to the end of next element. That is, in general, there are 
two springs attached to the two previous nodes. Fig. (3.2) shows the two-spring 
configuration model. 

Let Ksi be the stiffness of the spring attached at the node B and Ks 2 be the stiffness 
of the spring attached at the node C. 

S'l'e determined under the assumptions of Sec. 3.2 where each one of them 
acting alone and is enough to have crack up to node Ci. As can be seen, either Ksi 
or Ks 2 is the value of the spring stiffness attached at B or C necessary for the crack 
propagation upto Ci- That is any one of the springs will be able to sustain the crack 
upto Cl. 

In the present model, to simulate a smooth propagation of crack without any sudden 
jumps in the stresses etc., two springs act at any instant in general (or atleast one 
spring when the crack is upto an element boimdary). Therefore, while modelling 
the crack propagation with two springs, it is necessary, their stiffness values have to 
be halved, i.e, K 51/2 and Ks 2/2 for the sake of further analysis. 

That is, if the crack is upto Ci (within the element CD) as shown in Fig. 3.2(b), Ksi 
will be the stiffness value of the spring attached at B. Based on the model developed 
in Sec. 3.2, Ksi can be given by, 

Ksi — ^ ^ (3.16) 

where, jFCsi is the stiffness of the node B of the plate in the y-direction, for b.c. s 
shown in Fig. 3.2(b), fi is the displacement ratiom of node B given by, 

/, = iiL (3.17) 

Uoi 

where ui is the displacement in the y-direction of node B, and uqi is the displace- 
ment of node B when the crack tip reaches the end of ne«t element, i.e, D. This 
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also corresponds to the situation when the spring stiffiiess, Ksi becomes zero. As 
explained in Sec. 3.2, fi can be related to the crack length as, 


fi — ai+ Ci(ai — 


(3.18) 


where 


= A±l. 

di + 


(3.19) 


and C\ is a constant estimated from quasi-static crack propagation behaviour, 
similarly, Ksi will be the stiffness value of the spring attached at node C. Ksi is 
given by, 

Ks2 = K,i (3.20) 

where J\o 2 is the stiffness of node C of the plate in y-direction for the b.c.’s shown 
in Fig. 3.2(d), fi is the displacement ration of node C, given by 


/2 = — (3.21) 

^i02 

where ui is the displacement in the y-direction of.node C, and Uqi is the displacement 
of node C when the crack tip reaches the end of neat element, i.e, E, fi can be related 
to the crack length as, 

/2 = 02 + Ci{(ii — cif^) ( 3 . 22 ) 


where 


as shown in Fig. 3.2(c). 


02 = 


6, 

di + da 


(3.23) 


3.4 Closure ^ 

r 

/ 

Chapter 3 described in details of variation of various aspects of the proposed model. 
The formulation is done under th^ assu mption of one sp ring is actin g and it is 
generalised to the^p^e^e^odel. This model will be investigated with hv^pothetical 
data as well as experimental data in Chapter 4. 
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Figure 3.2; Crack opening scheme in present stiflEness release model 
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Chapter 4 

RESULTS AND DISCUSSION 


4.1 Introduction 

In this Chapter results of the dynamic crack propagation analysis are presented using 
the methods described in the previous Chapters. First the method is validated for 
both quasi-static case and dynamic case. Later results are also presented for the 
available experimental data for crack propagation. 

Various aspects of the results involve, 

• Determination of C for quasi-static crack propagation. 

• Crack propagation at constant speed under dynamic input force pulse. 

• Crack propagation at constant speed under time dependent input displacement 
(based on experimental data). 

4.2 static crack problem 

The F.E. code is first validated for a crack problem under static case. For this, 
a plate with center crack subjected to uniform tension over edges parallel to the 
crack axis as shown Figure (4.1), is solved for static case. Due to symmetry only 
one-quarter of plate is taken into account as shown in Fig. (4.2). 

The material properties and plate geometry are as follows; 
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• Modulus of Elasticity, E = 210 Gpa 

• Poisson’s ratio, v — 0.3 

• Density, p = 7840 Kgjm^ 

• Length of the plate, L - 0.104 m 

• Width of the plate, W = 0.04 m 

• Thickness of the plate, B = 0.024 m 

• Crack Length, 2a = 0.024 m 

• Virtual crack extension, Aa = 0.001 m 

• Far field stress, a = 

Mesh details are as follows: 

• Size of element in X-dir =1.0 mm 

• Size of element in Y-dir = 1.0 mm 

Theoretical value of energy release rate 4.084 J/m^ 

Energy release rate using present analysis 3.9426 Jfm^ 

Thus, percentage error = 3.46% 

Thus, it can be observed that present analysis is in good agreement with theoretical 
values (with in 3.5%). 

4.3 Dynamic Crack Propagation 

4.3.1 Determination of C 

p"' C found for a wedge loaded DCB specimen as shown in Fig. (4.3). Due 

to symmetry, only the upper half of the specimen is modelled by finite elements as 

shown in Fig. (4.4). 

Material and geometric properties are as follows. 
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• Modulus of Elasticity, E = 210 Gpa 

V 

• Poisson’s ratio, u — 0.3 

• Density, p ~ 7840 Kpjm^ 

• L(uigtli of the plate, L = 0.06 m 

• VVidtli of the plate, W = 0.005 m 

• Thickness of the plate, B = 0.024 m . 

• Crack Length, a = 0.03 m 

• Prescribed displacement, u = 0.025 mm 

The DCB siHicimo.n is analysed by using 4-noded isoparametric element. Details are 
as follows: 

• Number of elements = 600 

• Number of nodes = 671 

• Size of element in X-dir =1.0 mm 

• Size of element in Y-dir = 0.5 mm 

Energy release rate variation due to crack propagation is as shown in Fig. (4.5). 
It is observed that, for C = 0.0 though G is continuous within element, it is not 
continuous from element to element. It is found by trial and error for G = 0.62, G 
is continuous within element and also across element boundaries. 

This C value may not be suitable for dynamic case. However, this forms the initial 

guess value. 
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4.3.2 Effect of change in C value 

As obtained earlier for quasi-static case, value of constant C for smooth variation 
of energy release rate is (7 = 0.62. The same value when used for dynamic case, 
of velocity 1500 m/s the energy release rate curve is as shown if Fig. (4.6). Thus 
it can be seen from Fig. (4.6) that there are large oscillations at the inter element 
boundaries. It clearly shows that C value of the quasi-static case is not exactly 
suitable for the dynamic case and has to be adjusted. C value is adjusted so that 
the energy release rate curve becomes smoother. Effect of change in C value on 
energy release rate curve is shown in Fig. (4.7), for the crack velocity of 1500 m/s. 

It can be shown that for C = 0.4 give rise to less oscillations in G in comparison to 
other values of C and it is chosen that C = 0.4 for this case. Similarly for each case 
C value is found separately such that there will be less oscillations in the energy 
release rate curve. 

4.3.3 Hypothetical Problem 

For the analysis purpose a DCB specimen as mentioned in Sec 4.3.1 is considered. 
The problem is symmetric and only one half is considered for descritization and 
analysis. A short pulse is applied as shown Fig. (4.8). The pulse can be modeUed 
by F(1 - cos{ujt)) whose details are as follows; 

• Pulse duration P = 8 

• Force amplitude Fq = 50 N 

• Frequency of pulse / = 125 KHz 

Time step At depends on the element size and wave velocity. Crack propagation 
starts after the force pulse of Sij,s is completely applied. 

This problem is analysed for 3 different crack velocities and results are compared 
with that of Deshmukh [2000]. Same mesh as mentioned in Sec. (4.3.1) is used for 
all the cases. Details of time step and number of iterations are as follows: 
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Case 1; 


• Crack velocity = 1250 m/s 

• Crack initiation time = 8 ij,3 

• Time step, At = 0.05333 ij,s 
Crack is allowed to open up to 8 elements. 

Comparison of present model and Deshmukh [2000] is shown in Fig. (4.9). From 
the figure it is observed that the energy release rate curve of the present model is 
passing through the mean of that of Deshmukh [2000]. Oscillations are less in the 
present model. 

Case 2: 

• Crack velocity = 1500 m/s 

• Crack initiation time = 8 /rs 

• Time step, At = 0.06667 ns 

Crack is allowed -to propagate up to 8 elements. 

Energy release rate curve is shown in the Fig. (4.10). Oscillations in the present 
model are very less. 

Caise 3: 

• Crack velocity = 1750 m/ s 

• Crack initiation time = 8 /rs 

• Time step, At = 0.05714 

Energy release rate curve is shown in the Fig. (4.11). From the figure it is observed 

that oscillations are more in that of Deshmukh [2000]. 

^ — ^ ./ 

I 
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From the above 3 ca^es it is observed that present model even at high velocities also 
giving fairly smooth curve, but it is not the case of previous model. In the previous 
models when the crack propagates to the next element, the spring stiffness Ks is 
made zero and the plate is released suddenly, so there are oscillations at the element 
boundaries, but in the present model there is always atleast one spring is acting, so 
plate is not released suddenly at any time. 

4.3.4 Experimental data 

The present method is applied to determine the djmamic G for the experimental 
data as reported by Verma [1995] of a DCB. Material and geometric properties are 
as follows, 

• Modulus of Elasticity, E = 210 x 

• Poisson’s ratio, u = 0.3 

• Density, p = 7840 Kgjm? 

• Length of the plate, L = 0.16 m 

• Width of the plate, W = 0.0028 m 

• Thickness of the plate, B = 0.024 m 

Details of the mesh are as follows, 

• Number of elements = 4000 

• Number of nodes = 4411 

• Size of element in X-dir = 0.4 mm 

• Size of element in Y-dir = 0.28 mm 
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Crack length, crack initiation tune, velocity history and time history of 
displacement were found from experiment. This data is used as a input to the present 
code, and corresponding dynamic energy release rate, G is found. For four different 
cases, details are given in Table 4.1. Cantilever end deflection for experiment one is 
shown in the Fig. (4.12), similarly for experiment two, for experiment three and for 
experiment four it is shown in Fig. (4.14), Fig. (4.16) and Fig. (4.18) respectively. 
For this data the results are shown in Fig. (4.13), Fig. (4.15), Fig (4.17) and Fig. 
(4.19). Figures also shows the comparision with that of Siva Reddy [1997] 
Initiation toughness and propagation toughness values axe shown in the Table 4.2 
for different experiments. 

It can be seen from the figures that energy release rate is decreasing gradu- 
ally and fluctuations axe much less compared to the results reported in Siva Reddy 
[1997]. Propagation toughness for Elxpt No. 1 is 40% of initiation toughness, for 
Expt No. 2 is 30% of initiation toughness, for Expt No. 3 is 20% of initiation 
toughness and for Expt No. 4 it is approximately 1.5% of initiation toughness. It 
shows that at higher velocities initiation toughness is decreasing and propagation 
toughness is less percentage than that of initiation toughness. This shows at higher 
velocities the propagation toughness is less. 
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Expt No 

Crack length (mm) 

Crack initiation time {fj.sec) 

Crack Velocity {m/sec) 

1 

42 

44 

1200 

2 

42 

45 

1350 

3 

41 

43.5 

1500 

4 

42 

43.5 

1750 


Table 4.1: Details of experimenital data 


Expt No 

Gini 

Gprop 

1 

53 

20 

2 

60 

18 

■ 3 

50 

10 

4 

48 

7 


Table 4.2; Initiation and Propagation toughness for experimental data. 
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Figure 4.2: Plate with center crack (Half) 
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Figure 4.3: DCB Specimen (Full) 
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Figure 4.4; DCB Specimen (Half) 
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Figure 4.5: Energy releese rate vnriation for quasi-static crack propagation 
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Energy release rate (J/nr) 



Figure 4.9: Energy release rate variation for crack velocity of 1250 m/s 

(Hypothetical) 
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Figure 4.10: Energy release rate variation for crack velocity of 1500 m/s 

(Hypothetical) 
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Figure 4.11: Energy release rate variation for crack velocity of 1750 m/s 

(Hypothetical) 
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Figure 4.12; Displacement of cantilever end (Expt 1) 
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Energy release rate (J/nr) 



Figure 4.13: Energy release rate variation for case 1 of experimental data 
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Cantilever End Deflection, mm 



Figure 4.14: Displacement of cantilever end (Expt 2) 
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Cantilever End Deflection, 



Figure 4.16: Displacement of cantilever end (Expt 3) 
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Figure 4.18: Displacement of cantilever end (Expt 4) 
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Chapter 5 

CONCLUSIONS AND SCOPE 
FOR FUTURE WORK 

5.1 Conclusions 

A new model is proposed to investigate high speed crack propagation. This model 
involves further modification of stiffness release model and based on the results and 
discussion in Chap 4, the following conclusions can be drawn. 

, 1. Model gives fairly smooth and more stable variation of energy release rate. 

2. Experimental results reveal that unlike force release model energy release rate 
doesn’t drop to a low value in very few steps but falls gradually which is a 
more practical result. 

5.2 Scope for Future Work 

1. Present model can be further modified to get still better results. 

2. Present method can be extended for Mode II and mixed mode crack propaga- 
tion. 

3. This method can also be extended for 3-Dimensional crack propagation prob- 
lems. 
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